This invention relates to the field of geophysical prospecting and, more particularly, to a method for determining travel times and travel time maps for migrating seismic data.
In the oil and gas industry, geophysical prospecting techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a seismic energy source is used to generate a seismic signal which propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflections are recorded by seismic detectors located at or near the surface of the earth, in a body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the location of the subsurface reflectors and the physical properties of the subsurface formations.
In one form of geophysical exploration for natural resources, a seismic wavefield is radiated from a source point at or near the surface of the earth. Initially propagating as a spherically expanding wavefront, the radiation insonifies the various earth layers which usually offer an acoustic impedance mismatch at the layer boundaries due to variations in rock density and acoustic velocity. The wavefield is reflected from the respective layer boundaries to return to the surface where the mechanical earth motions due to the reflected wavefield are converted to electrical signals by transducers. The signals which comprise the seismic data, are stored on archival storage media for future processing.
It is the object of seismic studies to produce an image of a volume of subsurface earth formations in a region of geological and/or economic interest. For an isotropic horizontal stratum with constant velocity, the elapsed time between wavefield emission and wavefield reception at a receiver near the source, multiplied by half the average velocity, is the depth of the incident point of the reflected wavefield on the stratum that lies directly beneath the midpoint between the source and receiver.
If a reflector is tilted or the velocity is spatially variable, that simple relationship no longer holds; the incident point is shifted laterally up-dip relative to the source/receiver midpoint. Proper mapping of dipping or tilted reflectors, requires migration of the wavefields originating from those dipping strata.
Complex subsurface conditions preclude simple stacking of seismic records. Complex subsurface conditions scatter seismic waves in many directions. Simple stacking may position subsurface features in the wrong locations. Migration of seismic records can provide more accurate location of subsurface features. Migration involves geometric repositioning of return signals to show an event (layer boundary or other structure) in its proper location. A travel time map is an integral part of prestack migration processing. One well-known migration technique is Kirchhoff depth migration.
For prestack depth imaging, Kirchhoff migration is often chosen for its low computational cost relative to other methods. Kirchhoff migration requires use of wavefront travel time generators of any one of several well-known types that are based on ray tracing. Ray-tracing methods are useful in complex geological structures. Ray-traced travel times permit use of both first-arrival data as well as maximum-energy arrivals, which latter data produce superior imagery. However those methods are very slow and greedy of computer processing time.
The computation of wavefield travel times at selected output points is a key element in successful Kirchhoff depth migration of seismic data. For complex geology, wavefield travel times to the output points are generated by ray tracing from the source. Ray-traced travel-time generation methods produce travel times for all wavefront arrivals at an output point.
One method is taught by U.S. Pat. No. 5,229,938, issued Jul. 20, 1993 to Shein-Shen Wang et al. Here is taught a method for obtaining two-way travel times for source and receiver pairs that includes the steps of determining a set of one way travel times for each source to a plurality of image points and a set of one way travel times for each receiver to a plurality of image points. Ray sets are generated for both sources and receivers. Travel times from a source position.to image points are computed by two-point interpolation using the ray sets. Two-way travel time is computed by summing two sets, one set each for the source and receiver positions. A two-way travel time set is obtained for a particular source and receiver combination for all imaging points.
Eikonal (finite difference) travel time generators are very fast and do not produce shadow zones. Finite difference travel time generators will always pick the first arrival travel times whereas with ray-traced travel-time generators, the desired portion of the wave front must specifically be selected.
A well-known finite-difference travel time generator is disclosed in a paper published in the Bulletin of the Seismological Society of America, v. 78, n. 6, December 1988, pp 2062-2076, by John Vidale. Here the travel times of the first arriving seismic waves through any velocity structure can be rapidly computed on a multi-dimensional grid by finite-difference point-to-point extrapolation. Wavefronts are tracked instead of rays. Refracted waves are properly treated and shadow zones are filled with the appropriate wavefront segments. This scheme is very fast and is useful in tomographic inversion and Kirchhoff migration in geologic section characterized by smooth lateral velocity gradients.
Travel times, from which a travel time map or travel time table may be made, can be directly computed by ray tracing through a specified velocity model. A velocity model may be provided by well-known methods. A velocity model is generally considered to be a representation of the seismic wavefield velocity structure of the earth. The velocity model is a 2D or 3D array describing the distribution of velocities on a grid within the area or volume of interest. The grid may be a Cartesian (x-y-z) grid, although spherical, polyhedral or other grids may also be used. Determining a suitable velocity model for a given seismic data volume or earth structure is well known in the art and will not be discussed here in detail.
A bundle of rays emerging from a source location at the surface (of an earth velocity model) can be propagated down into the earth velocity model and traced through the subsurface while accounting for ray bending caused by changes in velocity gradient and refraction at layer boundaries with velocity contrast. Reflections points along each of the raypaths are identified as the intersection points of the rays with the layer boundaries. The travel time from the source location at the surface to a reflection point in the subsurface is then calculated by integrating the elements of distance along the raypath divided by the velocity associated with that element. By applying reciprocity, the travel time from a receiver location at the surface to a reflection point in the subsurface can be computed in the same manner. Finally, for a given source-receiver pair at the surface and a reflection point in the subsurface, the total travel time is computed by adding the travel time from the source to the reflection point the travel time from the reflection point to the receiver.
FIG. 1A is an example of a velocity model with one set of 4 raypath segments for a simple 3-layer earth model. Each of the 3 layers, V1, V2 and V3 have velocities such that the velocity of layer V1 is different from V2 and V2 is different from V3. A source location S at the earth surface from which seismic energy as a ray segment (represented by an arrow) propagates through the upper model layer V1. The next ray segment then takes a different path as it traverses the next earth layer V2. The total travel time from the source S to the reflection point M is the sum of the first to ray segments. At the interface of the second layer V2 and third layer V3 the ray reflects at reflection point M back toward the surface with velocity V2 and enters the upper layer V1, traverses this upper layer V1 and then arrives at receiver location R. For this source-receiver configuration, the total source to receiver travel time is computed by adding the ray segment travel times along the 4 individual segment paths.
Computed travel times may be mapped to grids that have a node structure analogous to velocity models. These travel time maps are a 2D or 3D arrays representing the distribution of travel times on a grid within the region of interest. These grids of nodes are analogous also to velocity models in that they may be Cartesian (x-y-z) grids, spherical, polyhedral or other grids.
FIG. 1B is an example of a 2D computation grid that may be used to map travel times from velocity models. All the circles, including S and R, in FIG. 1B represent grid points GP that represent nodes for which a travel time may be evaluated. The grid extends laterally in both directions x as indicated by the arrows, as well as in depth as indicated by the arrow in the z direction.
Travel times may be computed from a finite-difference solution to the eikonal equation or by ray tracing as described previously. However, both methods compute the travel time for an asymptotically high-frequency wave. These high-frequency methods do not map any finite frequencies accurately because the waves will pass through spatially varying velocity regions in frequency-dependent ways.
Kirchhoff migration produces an image by summing recorded reflection seismic data at image pixels, extracting a data sample for the sum at the estimated time required for the wave to travel from source to image point to receiver. Such travel times have been estimated by some solution to the high-frequency eikonal equation for over a decade. The first reported use of an xe2x80x9cupwindxe2x80x9d finite-difference method to solve the eikonal equation was by Wiggins, et al, 1986. The problems that arise in using high-frequency travel times in Kirchhoff depth migration are well known and much discussed since the early 1990s (Van Trier and Symes, 1991; Vinje, et al, 1993; Geoltrain and Brac, 1993; Audebert, et al, 1994a). Biondi, 1992, tried to improve on the high-frequency methods by performing a frequency-dependent smoothing of the velocity model.
Amplitude preservation during Kirchhoff migration is computationally expensive. For amplitude preservation, it is necessary to weight the migrated traces with suitable migration amplitude factors. Reference may be made to the article by Tygel, published in Geophysics, vol. 59, No. 12 of December 1993, which is a survey of the techniques employed. For example, it may be necessary to calculate a number of weights equal to the number of samples involved in the migration, and this is computationally expensive.
In a 1996 dissertation at Stanford, Dave Nichols computed monochromatic 2D Green""s functions and performed 2D migrations in an attempt to avoid limitations of the high-frequency assumption. A Green""s function for a given differential equation is the solution to the inhomogeneous equation with a spatial delta function as the source. Nichols computed monochromatic 2D Green""s functions using a finite-difference method in polar coordinates. Nichols was looking for a travel time to compute an approximation to the Green""s function that was a delta function in time. By treating a monochromatic Green""s function as approximately a sine wave, multiple frequencies may be computed and the resulting sine waves superimposed together as a peaked function. The function peak may be chosen as the location of the delta function.
The prior art discloses high-frequency methods and systems that do not accurately compute travel times through earth models with varying velocity structure. The prior art discloses methods and systems to compute migration amplitude factors that are computationally expensive. Accordingly, there is a need for a method and system to compute travel times more efficiently and accurately than prior art methods. There is a need for a system and method to compute migration amplitude factors efficiently to accurately migrate seismic images and represent earth structure.
A method for determination of seismic wave travel times from a surface location to at least one subsurface computation point (CP) within a plurality of subsurface location points. A velocity model of a region is defined. A monochromatic wavefield is computed, using the velocity model, for the at least one CP. An unwrapped phase of the seismic wave is determined at a neighboring computation point (NCP) proximate to at least one subsurface CP. The travel time of the seismic wave to the at least one CP is determined using the unwrapped NCP phase and the computed monochromatic wavefield for the at least one CP. Travel time maps and table may be created from the determined travel time values.